#!/usr/bin/perl
#### Genome_mutator.pl ###
## hart.steven@mayo.edu
use Getopt::Long;
use Benchmark;

my $inputChr = '';
my $outputVCF = 'snp_mutator_output.vcf';
my $mutations = 50000;
GetOptions(
	'f=s' => \$inputChr,
	'o:s' => \$outputVCF,
	'm:i' => \$mutations);
#print "intake options INPUT:$inputChr \nOUT : $outputVCF \nNUMMUT :$mutations\n";
if($inputChr eq '')
{	
	die "No Chr File Provided! \n"
}
print "Output: $outputVCF\n\n";
sub usage {
	print "usage: SNP_mutator.pl -f <Chr##.fa> \n";
	print "-o\tOutput File. Default = snp_mutator_output.vcf\n";
	print "-m\tTotal number of Desired Mutations\n";
	}
	
my @bases =("A","C","G","T");
our $bases_per_line = 50;
open (OUT, ">$outputVCF") || die "Can't write to output to $outputVCF!\n";
	
###################################
#### Start by Introducing SNVs ####
###################################
$insertions=int(($mutations-0.1)/10)+1;
$deletions=int(($insertions-0.1)/2)+1;
$random_places=($mutations+$insertions+$deletions);
#print " test: $insertions\n$deletions\n$random_places\n";
#die;
## If for whatever reason, numbers are used in place of ChrX,Y,M --> replaced
$inputChr=~s/chr23(.*)chr24(.*)chr25/chrX${1}chrY${1}chrM/;
## Find where the Genetic Code begins, avoid N's
$minimum_lines_to_skip =0;
$t0 = Benchmark->new;  ### ---------> Benchmark
open(INPUT,"<$inputChr") || die "Unable to Open $inputChr\n";
while(<INPUT>)
{
	next if($_ =~ /^\>/); ##Skip headers
	if($_ =~ m/^N/){$minimum_lines_to_skip++;}
	else{last;}
}
#while(<INPUT>){}
#my $max_lines = $.; ##Conitue to Last Line
my $max_lines =`wc -l $inputChr`;
$max_lines =~ s/\D//g;
#die "test $minimum_lines_to_skip $max_lines\n";
#print "Working Lines: $minimum_lines_to_skip - $max_lines\n";
$t1 = Benchmark->new;
$td = timediff($t1, $t0);
#print "Read Chr File took: ",timestr($td),"\n";
### Set initial Parameters #######################################
my $max_base = ($max_lines * $bases_per_line)-($minimum_lines_to_skip*$bases_per_line);	#find out how large the chromosome is
my $min_base = ($minimum_lines_to_skip*$bases_per_line);
#print "test1 min&max base :$min_base $max_base\n";
my @random_numbers=0;
my @random_numbers = map {random_int_between($min_base,$max_base)}(1..$random_places);#Choose a series of positions to mutate
print "test0 $min_base $max_base\n";
my @random_snps=@random_numbers[0..$mutations-1];
my @random_insertions=@random_numbers[$mutations..$mutations+$insertions-1];
my @random_deletions=@random_numbers[$mutations+$insertions..$random_places-1];

  @random_snps=sort {$a <=> $b} @random_snps;
  @random_insertions=sort {$a <=> $b} @random_insertions;
  @random_deletions=sort {$a <=> $b} @random_deletions;
print "test2 @random_snps \n@random_insertions\n@random_deletions\n";
#################################################################

$t2 = Benchmark->new;
$td2 = timediff($t2, $t1);
print "Generate Random Arrays took: ",timestr($td2),"\n\n";

seek INPUT,0,0;  #Reset to file beginning without Open/Close.
$.=0; ## reset line counter.
$chrom = <INPUT>;
my $chrid;
if($chrom =~ /chr(.+)/){$chrid=$1;}


###################################
#### Create SNV positions #########
###################################
$complete=0; 
my $n=0;
$target_line = lineOfSNP( $random_snps[$n] );
$last_line = lineOfSNP( $random_snps[-1] );
print "$random_snps[$n] $random_snps[-1] $target_line $last_line\n";
while(<INPUT>){
	if($. < $target_line){next;}
	elsif($. >= $last_line){last;}
	#print $.."\n";
	my $target_position= $random_snps[$n]-($target_line*$bases_per_line);
	print "test each line: $. $target_position $random_snps[$n] $_\n";
	if( $target_position > $bases_per_line){$n++; $target_line = lineOfSNP( $random_snps[$n] ); next;} ## math error handle
	@nucleo = split('', $_);
	if( !defined($nucleo[$target_position]) ){$n++; $target_line = lineOfSNP( $random_snps[$n] ); next;} ##just in case
	if( $nucleo[$target_position] eq "N" ){$n++; $target_line = lineOfSNP( $random_snps[$n] ); next;} ### Skip if base is N
	if( $nucleo[$target_position] =~ /\n/ ){$n++; $target_line = lineOfSNP( $random_snps[$n] ); next;} ### Skip if base is return character
	print "suc $nucleo[$target_position]\n";
	$base=randNonSelfBase($nucleo[$target_position]);	
	#print OUT "$chrid\t$random_snps[$n]\t.\t$nucleo[$target_position]\t$base\t.\t.\t.\n";
	print OUT "$chrid\t".$..":$target_line-$target_position:$n:$i\t$random_snps[$n]\t.\t$nucleo[$target_position]\t$base\t.\t.\t.\n";
	$complete++; $n++; $target_line = lineOfSNP( $random_snps[$n] );
	
}
print"SNVs requested: ".@random_snps."\tcompleted: $complete\n";

$t3 = Benchmark->new;
$td3 = timediff($t3, $t2);
print "SNV positions took: ",timestr($td3),"\n\n";

###################################
#### Now add insertions ###########
###################################
$complete=0; my $n=0;
seek INPUT,1,0;  #Reset to file beginning without Open/Close.
$.=0; ## reset line counter.
$target_line = lineOfSNP( $random_insertions[$n] );
$last_line = lineOfSNP( $random_insertions[-1] );
while(<INPUT>){
	if($. < $target_line){next;}
	elsif($. >= $last_line){last;}

	my $target_position= $random_insertions[$n]-($target_line*$bases_per_line);
	if( $target_position > $bases_per_line){$n++; $target_line = lineOfSNP( $random_insertions[$n] ); next;} ## math error handle
	@nucleo = split('', $_);
	if( !defined($nucleo[$target_position]) ){$n++; $target_line = lineOfSNP( $random_insertions[$n] ); next;} ##just in case	
	if( $nucleo[$target_position] eq "N" ){$n++; $target_line = lineOfSNP( $random_insertions[$n] ); next;} ### Skip if base is N
	if( $nucleo[$target_position] =~ /\n/ ){$n++; $target_line = lineOfSNP( $random_insertions[$n] ); next;} ### Skip if base is return character
	$insertion=randInsertionArray($nucleo[$target_position]);
	print OUT "$chrid\t$random_insertions[$n]\t.\t$nucleo[$target_position]\t$insertion\t.\t.\t.\n";
	$complete++; $n++; $target_line = lineOfSNP( $random_insertions[$n] );
	
}
print"Inserts requested: ".@random_insertions."\tcompleted: $complete\n";

$t4 = Benchmark->new;
$td4 = timediff($t4, $t3);
print "Insertions took: ",timestr($td4),"\n\n";

###################################
#### Now add deletions ############
###################################
$complete=0;my $n=0;
seek INPUT,1,0;  #Reset to file beginning without Open/Close.
$.=0; ## reset line counter.
$target_line = lineOfSNP( $random_deletions[$n] );
$last_line = lineOfSNP( $random_deletions[-1] );
while(<INPUT>){
	if($. < $target_line){next;}
	elsif($. >= $last_line){last;}
	
	my $target_position= $random_deletions[$n]-($target_line*$bases_per_line);
	if( $target_position > $bases_per_line){$n++; $target_line = lineOfSNP( $random_deletions[$n] ); next;} ## math error handle
	@nucleo = split('', $_);
	if( !defined($nucleo[$target_position]) ){$n++; $target_line = lineOfSNP( $random_deletions[$n] ); next;} ##just in case	
	if( $nucleo[$target_position] eq "N" ){$n++; $target_line = lineOfSNP( $random_deletions[$n] ); next;} ### Skip if base is N
	if( $nucleo[$target_position] =~ /\n/ ){$n++; $target_line = lineOfSNP( $random_deletions[$n] ); next;} ### Skip if base is return character
	$deletionLn=randDeletionArray($target_position);
	$deletion  = substr $_, $target_position, $deletionLn;
	chomp($deletion);
	print OUT "$chrid\t$random_deletions[$n]\t.\t$nucleo[$target_position]\t$deletion\t.\t.\t.\n";
	$complete++; $n++; $target_line = lineOfSNP( $random_deletions[$n] );
	
}
print"Deletes requested: ".@random_deletions."\tcomplete: $complete \n";

$t5 = Benchmark->new;
$td5 = timediff($t5, $t4);
print "Create Deletions took: ",timestr($td5),"\n\n";

close OUT;

########################################################################################################################
### SUBROUTINES
########################################################################################################################
sub random_int_between {
		my($min, $max) = @_;
		# Assumes that the two arguments are integers themselves!
		return $min if $min == $max;
		($min, $max) = ($max, $min)  if  $min > $max;
		return $min + int rand(1 + $max - $min);
		}
		
sub lineOfSNP {
	my $value = $_[0]/$bases_per_line;
	my $target_line=int($value);
	return $target_line;
}
sub randNonSelfBase{
        @self = ( uc($_[0]) );
        do
        {
                $b = @bases[rand @bases]
        }while($self[0] eq $b);
        return $b;
}

#sub randNonSelfBase{
#	@self = ( uc($_[0]) );
#	my @remaining = grep{ not $_ ~~ @self } @bases;
#	$b = $remaining[rand @remaining];
#	return $b;
#}
sub randInsertionArray{
	my $max_base = 50;
	my $min_base = 1;
	@insertion_length1 = map {random_int_between($min_base,$max_base)}(1..2);
	my $max_base = 5;
	my $min_base = 1;
	@insertion_length2 = map {random_int_between($min_base,$max_base)}(1..60);	# Make small insertions 30x more likely than large ones
	@insertion_length3 =(1) x 120;
	@insertion_length=(@insertion_length1,@insertion_length2,@insertion_length3);
	my $insertion_length=$insertion_length[rand @insertion_length];
	my $insertion = join '',map $bases[rand @bases],0..$insertion_length;
	return $_[0]."".$insertion;
}

sub randDeletionArray{
	my $max_base = 50-$_[0];
	my $min_base = 1;
	@deletion_length1 = map {random_int_between($min_base,$max_base)}(1..5);
	@deletion_length3 =(1) x 5;
	@deletion_length=(@deletion_length1,@deletion_length3);
	my $deletion_length=$deletion_length[rand @deletion_length]+1;
	return 	$deletion_length;
}


